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ABSTRACT 

Curves in a family derived from powers of the polar coordinate formula for 
ellipses are found to provide good fits to bound orbits in a range of power-law 
potentials. This range includes the well-known 1/r (Keplerian) and logarithmic 
potentials. These approximate orbits, called p-ellipses, retain some of the basic 
geometric properties of ellipses. They satisfy and generalize Newton’s apsidal 
precession formula, which is one of the reasons for their surprising accuracy. Be¬ 
cause of their simplicity the p-ellipses make very useful tools for studying trends 
among power-law potentials, and especially the occurence of closed orbits. The 
occurence of closed or nearly closed orbits in different potentials highlights the 
possibility of period resonances between precessing, eccentric orbits and circular 
orbits, or between the precession period of multi-lobed closed orbits and satellite 
periods. These orbits and their resonances promise to help illuminate a number 
of problems in galaxy and accretion disk dynamics. 

Subject headings: stellar dynamics — galaxies: halos — galaxies: kinematics and 
dynamics 


1. Introduction 

In recent decades, the recognition that galaxies are surrounded by extended halos of 
dark matter has spurred a greatly increased interest in the structure of orbits in gravitational 
potentials that are shallower than the Keplerian 1/r potential (e.g., Binney & Tremaine 1987, 
Bertin 2000). Similarly, the mass distribution in galactic bulges and elliptical galaxies is also 
extended, and the potential shallower than the Keplerian one. The logarithmic potential, 


$ oc ln(r), is among the simplest and most studied examples of such a shallow potential, and 
will be used as a primary example in this work. The celestial mechanics of such potentials 
is not nearly as well developed as the centuries old literature of the classical Kepler-Newton 
potential, but there has been much recent progress (see Boccaletti & Pucacco 1996 and 
Contopoulos 2002). 

The orbits in most physically relevant potentials can be numerically integrated very 
quickly and accurately, except near singularities and resonances. Direct numerical integration 
has provided much information about orbit families in both general symmetric, and non- 
axisymmetric potentials relevant to barred galaxies (including non-axisymmetric logarithmic 
potentials, see e.g., Miralda-Escude & Schwarzschild 1989). 

There has also been much progress in the areas of analytic approximation and global 
properties of orbits in non-Keplerian potentials, and power-law potentials in particular. Con¬ 
cerning global properties, Stoica & Font (2003) have studied the topology of energy surfaces 
of the logarithmic potential and derived a complete qualitative description of orbits in the 
axisymmetric case. They also proved that orbits are centrophobic for the non-axisymmetric 
case as well. Valluri et al. (2005) have presented extensions of Newton’s apsidal precession 
theorem (Newton 1687) for power-law potentials, and in the process found no evidence for 
cases of zero precession at high eccentricities. The work below provides further support for 
that conjecture. 

With regard to analytic approximation, we note first of all that classical epicyclic ap¬ 
proximations are the most well developed (e.g., see Binney & Tremaine 1987, Dchnen 1999, 
and Contopoulos 2002). Beyond this, Contopoulos & Scimenis (1990) have used the multiple 
Fourier series method of Prendergast (1982) to find orbit families and to derive accurate orbit 
approximations for the symmetric and non-axisymmetric logarithmic potentials. Adams & 
Block (2005) have examined spirograph or epicycloid orbit approximations and found that 
they accurately fit a wide range of orbits in the Hernquist and other halo potentials (with 
only two “spirograph” circles). The spirograph or epicycloid approximation is formally dis¬ 
tinct from the epicyclic approximation, though similar in that they both can be extended to 
arbitrary accuracy. 

More recently, Touma & Tremaine (1997) developed a relatively simple symplectic map¬ 
ping for non-axisymmetric potentials that captures the behavior of orbits in these poten¬ 
tials. The map is based on the approximation that orbital evolution is driven by two pro¬ 
cesses: apsidal precession in the axisymmetric part of the potential, and torques by the 
non-axisymmetric potential, concentrated at apoapse. Because this symplectic map is based 
on apsidal precession, the Valluri et al. paper and the results below on precession complement 
it. 
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With this recent work the knowledge base on orbits in galaxy potentials (and similarly 
shallow potentials in galaxy groups and clusters) has become quite extensive. However, the 
basis of much of classical celestial mechanics (and many specific approximation schemes) is 
the simplicity of Kepler’s elliptical orbit solution. The conceptual and analytical foundation 
of galactic dynamics and other fields that use power-law potentials is inevitably weaker with¬ 
out such simple orbital solutions. Unfortunately, mechanics textbooks teach us that analytic 
orbital solutions only exist in a few special cases among the power-law forces (Goldstein 
1980, Danby 1988). The simple epicyclic approximation provides a versatile computational 
tool, but with limitations as a conceptual tool if multiple or large epicycles are needed for 
accurate calculations (e.g., Contopoulos 2002, section 3.1). The references above show that 
other accurate approximations have been developed, but they are complex. 

The purpose of this paper is to present a family analytic orbit approximations for the 
logarithmic and other spherically symmetric, or axisymmetric potentials. These orbit func¬ 
tions are simply powers of precessing ellipses, and so, very simple indeed. Nonetheless, for 
small-to-moderate eccentricities, they are also surprisingly accurate. For the logarithmic 
potential, the approximate orbits librate around accurate numerical solutions, but do not 
diverge from them over time. Equally surprising, these approximate solutions extend contin¬ 
uously across all members of a family of potentials that range from the Keplerian potential 
to potentials with slightly rising rotation curves. Although they do not provide complete 
analytic solutions, their simplicity and accuracy allow them to fill a large part of the role 
of simple analytic solutions in providing a readily understandable conceptual picture of the 
nature of orbits in these potentials and the relations between them. 


2. Equations for Core Plus Power-Law Halo Potentials 

2.1. Potentials 

In this section I present the basic equations and dimensionless parameters that will be 
used in the following sections. All of the potentials used in this work give accelerations 
described by the equation, 


. . -GM*r 
9 ^ “ (e 2 + r 2 ) n ’ 

3 

with M(r) = r M«, 
w (e 2 + r 2 ) n 

so M* = 2 n M(e)e 2n “ 3 , 


( 1 ) 
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where e is a core scale length and M* is the constant scale mass of the potential. The 
last equation relates the scale mass to the mass contained within the scale length. The 
acceleration of equation (1) has many useful limiting forms. First of all, when r e, the 
acceleration increases linearly with radius, as does the circular orbit velocity. Thus, we have 
a solid-body rotation at small radii. In the special case of n — 3/2 we have the well-known 
Plummer model. The potential corresponding to this acceleration is also related to those 
studied in general scattering theory (see e.g., Newton 1966). 

Secondly, when r 3> e, we have the power-law limit, g oc l/r n_1 . It is convenient to 
define n — 1 + 5, so that the exponent of r is 5 in this limit. Three interesting special 
cases are: 1. n — 0, 5 — —1, circular velocity v cir oc r, 2. n — 1, 5 = 0, circular velocity 
v cir = constant, and 3. n = 3/2, 5 = 1/2, circular velocity v cir oc r -1 / 2 . Clearly, these cases 
correspond to linearly rising, flat and Keplerian rotation curves, respectively. 

The potential and density equations corresponding to the acceleration of equation (1) 
are given by, 


P 




- GM » 

25 (e 2 + r 2 )* 5 ’ 


3M* e 2 + (±4p) r 2 
ItT (e 2 + r 2 ) 2+<5 


( 2 ) 


The latter equation shows that density is approximately constant within the core. This 
equation yields negative densities when 5 > 1/2, and r 2 > (3/(25 — l))e 2 , i.e., at large radii 
in potentials where the outer falloff is more rapid than the Keplerian one. Such potentials 
are not very relevant to galaxy dynamics. However, some examples of them are used below, 
so it is worth noting that the density is positive out to quite large radii for potentials where 
5 is only slightly larger than 1/2. 

Currently, both observations (e.g., Diaz et al. 2005) and numerical structure formation 
models support the idea that halo profiles are universal. Among the most popular forms are 
the Navarro, Frenk, & White (1997, henceforth NFW) and Hernquist (1990) profiles. The 
density functions of both these profiles have a charateristic scale length, but they are also 
singular at the center, so their scale length does not correspond to a core size as in equation 
(2). Power-law potentials that are singular at the center will also be considered below. 
However, there is evidence that the core density profiles of some galaxies are flat rather than 
singular (e.g., Donato, Gentile, & Salucci 2004 and references therein). Henriksen (2005) has 
recently argued that the NFW state is a partially relaxed state, that in isolation ultimately 
relaxes to a core profile with solid body rotation. (For other recent views on the origin of 
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halo profiles see the recent works of An & Evans 2005, Dehnen & McLaughlin 2005, and Lu 
et al. 2005.) 

At large radii the NFW and Hernqnist density profiles fall off as 1/r 3 . This is also true 
of the profile of equation (2) when <5 = 1/2. 1 would emphasize, however, that the goal of 
this work is not to study how well a certain class of orbital approximations applies to specific 
halo potentials, Rather it is to describe how this family of approximations works in a range 
of simple schematic potentials, such as those of equation (2). 


2.2. Dimensionless Variables 

It will be convenient to work in dimensionless variables. The dimensionless radius is 
x = r/e, and the dimensionless time is r = i/ r //- T ff is a characteristic dynamical time, i.e., 
the free fall time at r = e. In terms of these variables, the Newtonian equation of motion 
can be written, 


d 2 x —fix 
dr 2 (l + a; 2 )”’ 


( 3 ) 


with the parameters, 


„ _ GM *Tff _ V <t> r Tff 
I e 2 n * 1 e 2 

where (f) is the azimuthal coordinate, v<p is the azimuthal velocity, and h is the dimensionless 
specific angular momentum. 

To further simplify the equation of motion we use the usual Kepler case substitutions 
of the inverse radius, u = 1/x, for the radius and the azimuthal advance, d(p = ( dcj)/dt)dt for 
the time change. The resulting equation of motion is, 

d 2 u f u 2 \ 5 1 / u \ GM * 

d^ + u ~ c \TT^ 2 ) V(i + « 2 )v = ’ C = ~h^ 

Henceforth, we will abbreviate du/d<p as u'. All of the results below apply to limiting or 
special cases of this equation. 

Before proceeding to these results we should briefly discuss the magnitude of the constant 
c, which is the only parameter in equation (5). Using the definitions of the parameters /i 
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and h (eq. 4) we have, 


c = 


( 


GAL 


l £ 2n-4 


(rvj 


( 6 ) 


In the case of circular orbits this reduces to, 


c = 


(! + xliX 


xz 


(7) 


where x cir is the dimensionless orbit radius. We see that for x cir « l,c oc \/x\ ir , while for 
x C ir » l,c oc l/x^~ 2n , which decreases to zero at large radii for potentials with n < 1. At 
x C i r = 1, c = 2 n , and thus, for orbits of this size c is of order a few in most of the potentials 
of interest. 


3. The Logarithmic Potential 

3.1. p-ellipse Solution to First Order in the Eccentricity 

We begin with the case of the logarithmic potential. The equation of motion of the 
softened logarithmic potential is given by equation (5), with 5 = 0. In the case of large 
orbits, or negligible core softening length, the equation of motion reduces to 

uu" = c — u 2 . (8) 

Henceforth, in this and all sections except 3.3, where the softening is treated explicitly, we 
will take e = 1, x — r, and u = 1/r, as is the convention for power-law potentials. 

Orbits derived from numerically integrating this equation look like precessing ellipses. 
Like elliptical orbits they are bounded by two turning-points (see e.g., Binney & Tremaine 
1987, Valluri et al. 2005). These facts motivate the trial of a solution of the form, 

u = -/(</>) [! + ecos((l - b)<j >)], (9) 

P 

where / is an unknown function used to track the deviations from the ellipse which the re¬ 
maining factors describe. That is, e is the eccentricity, the b factor determines the precession 
rate, and the semi-major axis of the ellipse a is given by p = a(l — e 2 ). p is the semilatus 
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rectum (see e.g., Ch. 2 of Murray & Dermott 1999). (Valluri et al.’s summary of Clairaut’s 
rotating ellipse and other historical approaches to the lunar theory are also relevant.) 

The solution of equation (9) solves equation (8) up to terms of order e 2 , if we take / 
equal to the inverse square root of the ellipse. Specifically, this approximate solution is, 

1 1 

u — - [1 + e cos ((1 — b)<p) l 1 / 2 , with. c\ = —, and, b\ = 1 — V2. (10) 

p p 2 

(For details see Appendix A.) 

The second equality above gives the first order approximation c\ for the constant c 
of equation (8) in terms of p, or the semi-major axis a. However, p is now the “semilatus 
rectum” for an oval that equals the square root of an ellipse. We will call a curve like equation 
(10), that is a power of an ellipse, a p-cllipse, for (precessing) power ellipse. Note that the 
parameter e functions in some ways like the ellipse eccentricity, but not in all ways, see 
discussions following equations (11) and (23) and in the appendices. For convenience we will 
continue to call it the eccentricity (see Adams & Block 2005 for a discussion of eccentricity 
definition). 

Many of the geometric properties of ellipses are retained by p-cllipses. For example, the 
major axis of a p-cllipse from equation (10) is given by a simple function of p and e, 


2 a = r(0 = 0) + r(( 1 — b)(j) 


n) =p 


\/l + e + yj\ — e 


( 11 ) 


Then the combination CiO 2 is a simple function of e whose value only varies from 1.0-1.27 as 
e is increased from 0 to 0.7. 

The two elliptical focii are still important for p-ellipses. On the other hand, the relation 
between the true anomaly and eccentric anomaly angles, as conventionally defined in celestial 
mechanics (see Murray & Dermott 1999), is much more complicated for p-ellipses than for 
ellipses. 

As conic sections, ellipses are tilted planar cuts of circular cones. Similarly, p-ellipses 
can apparently be derived as tilted planar cuts of circular paraboloids given by the formula 
z 2 = ( x 2 -( -y 2 y^ 2+5 , though this is difficult to prove in general. In the case of the logarithmic 
potential, the resulting Cartesian equation for the <5 = 0 p-cllipse is a fourth order polynomial. 
The Cartesian equation for the p-cllipses is not the same as that for super ellipses or the 
Lame curves (defined by the equation ( x/a) n + ( y/b) n = c), nor the elliptic or hyper-elliptic 
curves (which are quadratic in the variable y). These families of curves have quite different 
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The approximate solution given by equation (10) is technically only good to first order 
in e, so we might expect it to have quite limited accuracy, except for very small values of 
e. However, Figures 1 and 2 compare the fit of equation (9) to a full numerical solution of 
equation (8) in the case of a moderate eccentricity of e = 0.3. Figure 1 shows that the orbital 
fit is excellent over several orbital periods. Figure 2 shows that the approximation does not 
systematically drift away from the numerical solution after many cycles, and in fact is very 
close except near periapse. This shows that the approximate value of the precession factor 
6 is surprisingly accurate. (Part of the reason for this accuracy is the slow variation of the 
constant c with eccentricity discussed above. These points will be discussed further below.) 


3.2. Second Order Approximation (p-ellipse plus epicycle) 

A clue to the lack of drift between the approximation (eq. (10)) and the solution noted 
in the last subsection is that the only second order terms that appear are proportional to 
(ce 2 ) sin 2 ((l—6)0) and (ce 2 ) cos 2 ((1—6)0). Moreover, there are no third or higher order terms. 
The fact that these are simple sine or cosine squared terms suggests that we could cancel 
them with more trigonometric functions in the approximate solution. Thus, we next consider 
approximations consisting of the sum of a p-cllipse and an elliptical epicycle (constant plus 
cosine term), 


u = — [ 1 + e cos ((1 — 6)0)] 1 ^ 2 + — [1 + e 2 cos ((1 — 6)0)1. (12) 

Pi P2 

As before, we substitute this into equation (8) and equate the sum of terms with common 
powers of ecos((l — 6)0) to zero. Additionally, we assume that e 2 oc e. We find that the 
resulting constant terms yield an expression for c, the first order terms yield the same value 
of 6 as in the previous subsection, and the second order terms yield a relation between e 2 /e 
and Pi!p\- If we further demand that equation (12) equals equation (10) at 0 = 0 (and 
that both equal the initial value w(0) of the full solution), then we obtain another relation 
between these two ratios and can solve for both. The results are, 



c 2 = 


, 6 2 = 1 - y/2, e 2 = Vl + e - 1 , and, 
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The subscript “2” is used on the constants b and c to indicate that these values are accurate 
to order e 2 . 

It is assumed that the added cosine term has the same precession factor as the p-ellipse, 
so it does not drift from the p-ellipse approximation. A fortuitous cancellation leads to a 
zero second order correction to the precession factor b. We will see below that this is not 
generally true in other potentials. 

This second order solution is plotted along with the numerical and first order solutions 
in the example of Figure 2, and we see that it doesn’t provide much improvement. The 
reason for this is that in going from the first to the second approximation, we have changed 
the higher order terms remaining in the orbit equation from two second order terms that 
partially cancel at most azimuths, to more than half a dozen third order terms, which do not 
generally cancel each other. Note that in this example: p\ = 1.93, P 2 = 1.73, and e 2 = 0.14. 

These third order terms have coefficients that are third order powers of sines and cosines, 
so they could presumably be eliminated by adding still more trigonometric terms to the 
approximate solution. We will not pursue the details, but it thus appears that we could 
improve the approximation with the addition of more epicycles or Fourier terms. Because 
the p-ellipse part of the solution is itself a good approximation, and because successive 
Fourier coefficients contain successively increasing powers of the p-eccentricity, the solution 
series should converge steadily for moderate eccentricities. 


3.3. The Softened Logarithmic Potential 

The accuracy of the simple p-ellipse approximate solution for orbits in the logarithmic 
potential is remarkable. In this and the following subsections we will see that it is not 
unique. On the contrary, it is useful over a range of interesting potentials. In this subsection 
we consider cases of equation (5) with 5 = 0, but with nonzero core length e (see eq. (1)), 
corresponding to the softened logarithmic potential. 

In this case equation (5) can be written as, 

uu" + u 2 — c = —u 2 (uu" + u 2 ), (14) 

where the left-hand-side is the same as equation (8), and the additional nonlinear terms have 
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been placed on the right-hand-side. (Note that in this subsection we return temporarily to 
the convention of measuring all length units relative to e, i.e, here u = e/r.) The term in 
parenthesis on the right-hand-side of equation (14) is also very similar to equation (8), which 
is a large part of the reason that the p-ellipse approximations continue to work well in this 
case. 

Upon substitution for u from equation (10), and equating zeroth and first order terms 
in ecos((l — b)<p) as above (see Appendix A), we derive new forms for the factors b and c 
appropriate to this potential, 



In the limit of orbit sizes much larger than the softening length (p 1), these expressions 
reduce to those of equation (10), as they must. Additionally, in this case, we have two other 
limiting cases. The first is the limit of small orbit sizes where the rotation curve is linearly 
rising. In that case, b\ ~ —1, which as we will shortly see is very interesting. The other 
special case is the exchange or rotation curve downturn region, defined by p ~ 1, where there 
are also interesting resonances. 

Before proceeding to discuss these cases it should pointed out that if we substitute 
the second order approximation of equation (12) into equation (14), we find that as in the 
previous subsection, there is no second order correction to the b factor. The c factor, however, 
does have a correction of order e 2 , 



Figure 3 illustrates the goodness of fit of the (first order) p-ellipse solution, for different 
size orbits. As we would expect from the relation b ~ — 1, panel a) shows that the analytic 
solution is a slowly precessing, nearly closed, oval when p is small (here p = 1/2). The panel 
also shows that the analytic approximation fits the numerical solution quite well. It librates 
around it, but returns periodically, just as in the case of the pure logarithmic potential. 
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The second panel of Figure 3 is much like Figure 1, which corresponds to the large p 
limit in the present case. That is, the solution is a more rapidly precessing, open p-ellipse, 
which the analytic approximation tracks quite well. It is interesting that we can see this 
behavior already at p — 1, the case shown. 

The third and fourth panels of Figure 3 show special cases of resonant orbits in the 
transition region. The third panel shows the case b = —2/3, corresponding to p ~ 1.2536. 
Both analytic and numerical curves close (or nearly close) after 3 times around the center. 
The fourth panel shows the case b = —3/5, corresponding to p ~ 1.6036, and we have more 
loops before closure. 

The factor b ranges over the interval (—1,1 — y/2). The infinite rational numbers in 
this interval correspond to closed (analytic) orbits. However, if we write a rational number 
m/n, there are very few in this interval with small values of nr or n, i.e., low order period 
resonances. Like the above examples most of these will have radii of p ~ 1. We reiterate 
that the value of b does not depend on the eccentricity up to at least second order. The 
examples of Figure 3 have a moderate eccentricity of e = 0.3. 

In sum, Figure 3 shows that the first order p-ellipse approximation provides a good fit to 
orbits in this potential, as in the case of the logarithmic potential. In the previous subsection 
we found that the second order approximation yielded only a small improvement. In the 
present case the second order derivation involves a very unwieldy, high order polynomial in 
Pi and P2 ■ Since this analysis produces no significant new results, I will not present it here 
(except for eq. (16)). 


4. Potentials with Rising or Falling Rotation Curves 
4.1. General p-ellipses 

In the preceding section we considered the special cases of bound orbits in the logarith¬ 
mic or softened logarithmic potentials, i.e., where 5 = 0 in equation (5). In this section we 
consider positive and negative values of 5, which correspond to falling and rising rotation 
curves. We will not, however, consider the most general form of equation (5). To simplify 
the analysis and the discussion, we will neglect softening and consider only pure power-law 
potentials. (We will also return to conventional variables, where u = 1/r.) 

In this case the equation of motion becomes, 


uu" = cu 26 — u 2 . 


( 17 ) 
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The simple p-ellipse solution in this case is, 


1 


u — - 


V 


[1 + ecos ((1 - b)4>)}* +s , 


(18) 


with the variables defined as above. Specifically, we see that when 5 = 0 we have the 
approximation of the previous section for the logarithmic potential, and when 5 = 1/2 we 
have the ellipse for the Keplerian case. Given that in both of these cases equation (18) is an 
excellent approximation, it is not completely surprising that the equation is also an excellent 
approximation for all intermediate values of 5. We will return to that point in a moment. 

First we follow the procedure of the previous section. We substitute equation (18) 
into equation (17) and eliminate all terms of zeroth and first order in e, by making the 
identifications, 


ci = ( -) , and, bi — 1 — \/2(l - 5). (19) 

Comparison to equation (10) shows great similarity, with the addition of some simple factors 
of 1 — 5. The quantities p, e, and 5 determine a p-ellipse, but it is also desireable to specify 
an orbit in terms of the conserved specific angular momentum and specific energy. Relations 
between these variables are presented in Appendix B. 

In Figure 4 we show some sample approximate and numerically integrated orbits with a 
range of 5 values, and with p = 1, e = 0.2 in all cases. A quick glance at this figure reveals a 
couple of basic points. First, the structure of comparable orbits varies a great deal over these 
different potentials. Second, the first order p-ellipse approximation is very good in all cases. 
If we look more closely at the individual panels we see that, as in the case of the logarithmic 
potential, the approximation librates a short distance from the numerical solution, only to 
return to it shortly thereafter. As we will discuss in the next section, this is generally true 
because the precession frequency implied by the second equality of (19) is exactly right for 
all the potentials (i.e., the same as in Newton’s theorem). 

We see that the orbit shown in Figure 4b is similar to that of Figure 1. The exponent 
5 in this case is half way between that of the logarithmic and Keplerian potentials, so the 
fact that the orbit is similar to the former, but more slowly precessing, reflects continuity 
between these solutions. Figure 4c is very close to the Keplerian case, and the precession 
is very slow indeed. If we cross to the other side of the Kepler case, as in Figure 4d with 
5 = 0.7, both approximate and numerical solutions appear more erratic. This example is in 
the range between the Kepler solution and the Cotes spiral solution at 8 — 1.0 in the present 
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notation. (The MATLAB routine ODE45 was used throughout this paper.) In any case, we 
don’t expect the p-cllipse approximation to be very useful for values of 5 of about or greater 
than 1.0. 

Figure 4a shows an example of going past the logarithmic potential to those with nega¬ 
tive values of 5. The exponent 1/2 + 5 of equation (18) becomes very small in the vicinity of 
5 = —1/2, so even with significant values of the eccentricity e, the radial excursion of orbits 
is very small. Nonetheless, the figure shows that the p-ellipse approximation continues to be 
good. 

The two panels of Figure 5 show sample orbits at still more negative values of 5. The 
first panel shows a cases where 5 = —1.1 is very close to that of the solid-body potential, and 
so is quite similar to Figure 3a. The second panel shows a still more extreme case (5 = —1.5), 
and here the deviations between approximate and numerical orbits become quite significant. 
Thus, it seems that the p-ellipse approximation is at best only qualitatively useful for values 
of 5 much more negative than the solid-body case. 

The normalized (by a factor of p 2 ) second order terms that remain when the p-cllipse 
is substituted into the equation of motion provide an estimate of the inaccuracy of that 
approximate solution. Technically, they provide an estimate of the inaccuracy of the u"/u 
value, but given the common periodicity of u and its derivatives, this is roughly equivalent to 
the relative inaccuracy in u. All the second order terms are proportional to (e cos((l — 5)</>)) 2 . 
If we take the average value of 1/2 for the cosine term, a first rough estimate of the error 
is e 2 /2, which, when compared to the numerical results, turns out to be about right in the 
vicinity of 5 = 0. 

The 5 dependence of this error estimate is given by the function: (1 + 25) (1 — 5)(5 3 + 
0.55 2 + 3.55 — 1). When 5 values vary from -0.5 to 0.5, the absolute values of the function are 
always less than 1.2, and mostly in the range 0.5-1.2, i.e., of order unity. For values 5 <C —0.5, 
the absolute value of the function rises rapidly (eventually as (1 + 25)5 4 ), explaining the 
diminished accuracy of the p-cllipse approximation. Within the domain of 5 values of the 
most interest (-0.5 to 0.5), eccentricity is the dominant factor in the error estimate. Section 
5.2 offers a way to derive more accurate p-ellipse approximations at large eccentricity. 

The approximate nature of the first order p-cllipses, and these limitations on their do¬ 
main of accurate application, distinguish them from a universal analytic solution. However, 
the latter evidently doesn’t exist, and the p-cllipses provide accurate orbit approximations 
in power-law potentials over most of the domain of interest for studies of galaxy dynamics. 
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4.2. Second Order Approximation for Power-Law Potentials 

As in the case of the logarithmic potential, we can attempt to eliminate terms of second 
order in e by adding the p 2 , e 2 terms of equation (12) to equation (18). In the more general 
cases of 5 ^ 0.0 there are more second order terms than in the logarithmic potential, so we 
might expect the correction to be more useful than in that case. 

However, direct comparison of the second order approximation to numerical integrations 
in cases like those shown in Figures 4 and 5 shows that the situation is not that simple. I find 
that for 5 values in the range 0 < 5 < 1/2 there is little difference between the p-cllipse and 
second order approximations, as in the case of the logarithmic potential. For —1/2 < 5 < 0, 
the second order approximation fits the numerical orbits significantly better than the first 
order p-ellipse. For <5 < —1/2, the amplitude of the second order approximation fits better, 
but its precession rate is apparently not accurate because the orbital phase drifts quite 
rapidly. We conclude that, as in the case of the logarithmic potential, there is not much 
reason to prefer the second order approximation in the case of general power-law potentials. 

The following equations define the second order solution, 


Pi 

P2 


bo = 1 


2(1-5) 


1 1/2 


1 + |(1- 45 2 ) 
\ l+<5 


_P2_\ e 2 


P1+P2 ) 


e 2 = (1 + e) 2 — 1, (initial condition ), 


( 1 + 2,5) (f-f) 

Hf ) 2 


1 ± (1+2,5) 


3 e 2 


+ 2(1-25) 


pi = I 1 + — ) c 2 ^- 4 ) , 
P2, 


( 20 ) 


where the variables e, c, and 5 are given, and the above equations are solved for pi,p 2 ,e 2 , 
and b 2 - The equation labelled initial condition is derived by equating the initial value of 
the first and second order approximations. 

The expression for b 2 in equation (20) has a slight dependence on e. This allows the 
second order approximation to drift relative to the first order solution. It also allows for 
the possibility of a rational value of the precession period relative to the orbital period, and 
closed orbits, at special values of e. We will consider this further in Sec. 5. 
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4.3. Potentials with Closed Orbits 

The second equality of (19) allows for values of 5 that yield rational values of b±, and thus, 
closed orbits for all low-to-moderate values of e in those potentials, in the first order p-cllipse 
approximation. These potentials are generally isolated, i.e., not dense in the mathematical 
sense, on the 5 axis, though there are exceptions. The existence of these closed orbits 
and their basic structure can also be deduced from the apsidal precession theorem and 
its generalizations. However, there is the question of how small the eccentricity must be 
for Newton’s theorem to be valid, and the generalizations of Valluri et al. (2005) to high 
eccentricities are complex. If p-ellipse formulae like those of eqs. (19) and (20) are sufficiently 
accurate, then they provide simple alternatives for guidance on these questions. We examine 
this question in this subsection. 

From the second equality of (19) we deduce that in order for b\ to have a rational value, 
we must have, 


5 


2 k 2 ’ 


( 21 ) 


where j and k are integers. Then, b\ = (k — j)/k. As a simple and generalizable example 
consider the sequence k — 1,2, 3,..., with j — k+ 1. The corresponding sequence of b\ values 
is b\ — — 1/k. Table 1 gives 6 values for some elements of this sequence. 

Figure 6 gives examples of orbits in this sequence with eccentricities of e = 0.2 and orbit 
sizes of p = 1.0. Specifically, the k = 2, 3, 4, and 20 orbits are shown. The first panel shows 
a very simple closed orbit, that with S = —1/8 is close to the logarithmic potential. The 
analytic and numerical solutions are very similar to each other. To keep the plot readable 
we only show a couple of times around the orbit. However, a plot like that of Figure 2 for 
many orbits shows no drift between the numerical and analytic curves. I will not attempt to 
prove that the orbits in this potential are closed, but it appears that for practical purposes 
they are. (This has been checked at other values of e as well.) 

The second panel of Figure 6 corresponds to a slightly falling instead of a rising rotation 
curve with 5 = 1/9. It is slightly more complex, but still closed, and with excellent agreement 
between the analytic and numeric curves. As before, there is no drift between the two curves 
after many orbits. 

The trends of more complexity in the closed curves, and good agreement between ana¬ 
lytic and numeric curves, continue in the third panel of Figure 6, where 6 = 7/32. 

The fourth panel, where 6 ~ 0.45, shows that things get very interesting as we approach 
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the Kepler case. The orbits become very slowly precessing p-ellipses that eventually close. 
As k advances to large values in this sequence, b { goes to zero, and according to equation (21) 
5 goes to 1/2. Thus, the sequence converges on the Kepler case. In contrast to the sparsity 
of potentials in this sequence in the region of the logarithmic potential, the convergence of 
this sequence shows that such potentials do densely populate the 5 axis at values just less 
than the Kepler value. These statements are also true for a second sequence defined by k — 
1, 2, 3,..., with j = k + 2, and so on. 

Thus, the properties of the p-cllipses, originally derived as simple analytic approxima¬ 
tions, have revealed families of potentials with (at least practically) closed orbits, like the 
Kepler potential, and unlike the logarithmic potential. Appropriately, these families cluster 
thickly on one side of the Kepler potential. Interestingly, these special potentials include 
the 5 = —1/8 and 1/9 potentials, but not the simple 5 = 1/8 and 1/4 potentials. This is 
explained by the precession law of the p-ellipses. 

I have avoided discussing the very simple case of k — 1, j — 2, and b\ — 8 — — 1. This is 
the first element of the sample sequence above, and is the solid-body potential. The analytic 
and numeric curves for this highly resonant case do not agree. This is true in spite the fact 
that we have shown good agreement in examples above, which have 5 values a little less than 
- 1 . 0 . 


5. Discussion and Applications 
5.1. Apsidal Precession and Integrals of the Motion 

It has been noted already that one of the secrets of the success of the first-order p-cllipse 
approximations is that their precession rates accurately match those of the true orbits in 
power-law potentials up to moderate eccentricity values. It remains to demonstrate explicitly 
that they satisfy Newton’s apsidal precession theorem, and consider large eccentricities. 

As to the former, we note that the minimum of u (apoapse) occurs when (1 — 6)0 = tt, 
the maximum (periapse) when (1 — 6)0 = 0, and the difference is, 


A 0 


7r 


( 22 ) 


which is the apsidal precession theorem for p-cllipses. Valluri et al. (2005) remind us that 
with an acceleration proportional to l/r N ~ 3 , Newton’s apsidal precession theorem gives 
A0 = n/y/N. When we compare this form of the acceleration to the one used in this paper 
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(with no softening), we find N = 2(1 — 5), so the two expressions for the apsidal precession 
are the same. This proves the that the first-order p-cllipses satisfy the apsidal precession 
theorem. (Note that with eq. (15) we also can readily derive the apsidal precession as a 
function of orbit size in the softened logarithmic potential.) 

Valluri et al. (2005) emphasized that apsidal precession depends on orbital eccentricity 
at large eccentricities, so the Newton formula and equation (22) are not correct in that limit. 
These authors derived integral expressions for the dependence of the precession rate on N 
and e for the power law potentials. The first equation of (20) provides a simple approximate 
formula for this dependence. 

Before comparing this simple formula to the Valluri et al. results, we should recall that 
the p-cllipse eccentricity we have used up to this point is not the same as the usual definition 
of the eccentricity (except in the case of the inverse square law). If an orbit is assumed to 
be approximately a precessing ellipse, then one definition of eccentricity is in terms of the 
periapse distance using the formula r mtn = a(l — e). We will call the eccentricity defined 
in this way the classical eccentricity, and denote it e a . The classical eccentricity and the 
p-ellipse eccentricity used above are related by the equation, 

(1 - e ) S+ ^ = 1 - e 0 . (23) 

The classical eccentricity is generally significantly less than the p-ellipse eccentricity, though 
they are nearly the same for potentials close to the Kepler potential. 

With that notational disparity accounted for, we return to the topic of the dependence 
of precession on eccentricity, and the comparison of Valluri et al.’s results to the 62 expression 
in equation (20). First of all, it should be emphasized that we have seen that the second- 
order approximation described above was not found to be highly accurate in many cases. 
Valluri et al.’s semi-analytic apsidal precession angles, illustrated in their Figure 2, are more 
accurate than those derived from equation ( 20 ). 

Nonetheless, the u — <fi plot in Figure 7 shows an example where the precession of a 
numerical orbit solution after 10-20 orbits is much better fit by the second order precession 
rate than the first order rate. (Their libration amplitudes are comparable.) Equation (20) 
does seem to be an improvement over the b\ formula of equation (19). The example shown 
in Figure 7 has e = 0.95 and e Q = 0.84, so it is quite eccentric. The example in Figure 7 has 
5 = 1/9, a closed orbit value, but at this large eccentricity the orbits are not closed (see Fig. 
8 ). 

The 6 2 expression of (20) also captures and helps us to understand many of the quali- 
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tative features of Valluri et al.’s results. First of all, the dependence on eccentricity is quite 
weak in equation (20), and so, it only becomes important for large values of e (unless 5 is 
very large). Indeed, either equation (20) or Valluri et al.’s graphs show how good an approx¬ 
imation the constant precession rate b\ is for given values of e and 5 (and thus overcome the 
primary objection to the use of the apsidal precession theorem up to modest eccentricities). 
Secondly, both formulations predict that the precession rate is constant at all eccentricities 
for the logarithmic potential. Thirdly, they both agree that the apsidal angle A0 increases 
(decreases) with eccentricity when 5 < 0 (6 > 0). In fact, equation (20) provides a fair 
quantitative approximation to the plots of Valluri et al. at values of 5 not far from zero. 

Valluri et al. also state that their work gave no evidence for “isolated cases of zero 
precession as e was increased.” The b 2 expression in (20) indicates that such cases might 
exist at values of 5 just slightly less than 1/2, but even if this prediction is confirmed, these 
cases would be very anomalous. Another somewhat arcane prediction of equation (20) and 
Valluri et al.’s results is that potentials very near the resonant ones discussed in the previous 
section should have a closed orbit at one specific, high value of the eccentricity. These 
potentials have smaller values of 6 than the resonant one for positive 5. 

The discussion of this subsection clearly relates to the classical integrals of the motion. 
The precession rate b (or the precession frequency divided by the orbital frequency), is an 
integral of the motion, much like the initial orbital phase angle which is often used as such 
(see Sec. 3.1 of Binney & Tremaine 1987). By indicating the presence of closed orbits, the 
expressions for b in (19) and (20) tell us (approximately) when these integrals are isolating. 


5.2. An Alternate Approach to p-ellipse Orbit Approximation 

5.2.1. Better Appi'oximations at Larger Eccentricities 

The p-cllipse orbit approximations above were derived from an analytic perturbation 
analysis, with some implicit constraints. In this section we learn more about these approxi¬ 
mations by relaxing some of these constraints. Figure 7 shows how, in the case of an orbit 
of high eccentricity (5 = 1/9 and e = 0.95), the p-ellipse approximations drift in period 
relative to the numerical solution. (This is not surprising given the Valluri et al. (2005) 
results on apsidal precession.) This motivates an attempt to improve the p-cllipse solution 
by changing the value of b to match the numerical result. In the case shown in Figure 7 this 
means changing from the analytic value of b 2 = —0.3027 to -0.3558. 

According to equation (19) or (20), a different b factor corresponds to a different 6 
value. In this case (using eq. (19)), we go from 5 = 1/9 = 0.1111... to 6' ~ 0.0809. Another 
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difficiency of the p-ellipse approximation shown in Figure 7 is the radial range. With the 
revised value of 6 we can also revise p and e to fit the numerical maximum and minimum 
radius. The revised p-cllipse approximation is given by the equation, 

- = 1.1105 [1 + 0.7026 cos(1.35580)] 0 ' 5809 . (24) 

r 

This provides a very good fit to the orbit as shown in Figure 8, where the numerical orbit is 
shown by the solid line, and the approximation by the dashed line. 

To recapitulate, we have treated the p-cllipse quantities as free parameters, first adjust¬ 
ing the precession rate to match the numerical value, and then the orbital size and eccen¬ 
tricity to match the radial range. The 6 value was also changed to maintain the relationship 
of equation (19). This procedure seems to yield an highly accurate orbit approximation for 
such loop type orbits, even at quite high eccentricity. 


5.2.2. p-ellipse Transformation and Perturbation Approximation 

The result of the previous subsection can be viewed as achieving a more accurate fit 
to a high eccentricity orbit than the first-order p-cllipse by using a smaller eccentricity 
p-ellipse orbit taken from a different (power-law) potential. This suggests that there is 
an approximate mapping between p-ellipses, representing an approximate symmetry of the 
power-law potentials. The nature of this relation is that the orbits in a given potential 
at high eccentricity correspond to orbits of lower eccentricity in a different potential, with 
the same precession rate. Of course, the first-order p-cllipse approximation fits the lower 
eccentricity orbit better. 

As an aside, note that this relation suggests that periodic orbits may be found at high 
eccentricity in a potential with generally aperiodic orbits, when the correspondence is to an 
individual orbit in a periodic potential. 

However, the assumption that the first order equation (19) holds at high eccentricity is 
likely to be a rather crude ansatz. Another approach is to return to the perturbation analysis 
of Sec. 4.1, but now allow the 5 factor of the p-ellipse approximate solution to be different 
than that of the equation of motion, call it 5. We will still fix the b factor to the numerical 
result, because even a small precessional drift will result in a large difference between the 
approximate and true orbit. 

Then, requiring the cancellation of zeroth and first order terms in ecos((l — b)(f>) yields 
the following relation, 
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2(1 -5) 

1 - |(1 -4<5 2 )e 2 


(25) 


(Note the similarity to the expression for 62 in eq. (20)) 

For a given potential exponent 5 and precession parameter b, this gives a relation between 
the exponent 6 of the corresponding potential, and the eccentricity of the p-ellipse in that 
potential. Requiring that the approximation also fit the maximum and minimum radii of the 
numerical orbit (which can be derived from the conserved energy and angular momentum, 
see Appendix B) gives three relations for the variables <5, e, and p. Second order terms in 
the perturbation expansion do not cancel, but since the b value is correct, this procedure is 
better than the direct second-order approximation. The results for the case shown in Figure 
8 are: 5 = —0.003, p = 0.87, and e = 0.77. The values are very similar to those of equation 
(24). This perturbation approach is less ad hoc than that used to derive equation (24), and 
thus, provides a more solid basis for approximating high eccentricity orbits. 


5.2.3. Epicycloid Comparison 

A specific comparison to the epicycloid approximation of Adams & Block (2005, see their 
eqs. (20) and (21)) is also shown in Figure 8 , as a dotted curve. This curve was constructed 
by requiring that it have the same maximum and minimum radii, and the same precession 
rate as the numerical solution. It can be seen from Figure 8 , and the corresponding u — (j) 
graph that it deviates from the numerical solution a bit more at intermediate radii than the 
specially fitted p-ellipse. 

Overall, the simple versions of both approximations seem to be about comparably ac¬ 
curate in the case shown. It seems likely that the epicycloids, like the p-cllipses, can provide 
good orbit approximations in the potentials studied above. Adams & Block (2005) have 
considered the properties of the epicycloid orbits in general potentials, so I will not examine 
this issue in detail in this paper. The advantages of the p-cllipses include their simple form 
in polar coordinates, and the fact that their orbit parameters change in very simple ways as 
functions of the conserved quantities and the parameter 5 in power-law potentials. 
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5.3. Galaxy Collisions, Bars and Precession Period Resonance 

In this paper, I will not present any detailed applications of the p-cllipse theory devel¬ 
oped above. Indeed, since the orbits have long been studied numerically in many relevant 
cases we do not expect new applications, as much as new perspectives on old applications. 
Nonetheless, in this section I will highlight a couple of applications that seem worthy of 
further development. 

The first is simply the application of p-cllipses to obtain the general characteristics of 
satellite orbits in potentials with either rising or falling rotation curves. The simplicity of p- 
cllipse allows one to quickly produce a wide range of possible satellite orbits in any particular 
situation. In studies of interacting galaxies, one can then use the impulse approximation to 
estimate the effect of a close encounter on initially circular stellar orbits in the disk of the 
primary. The initial effects of such encounters are kinematic, with self-gravitational effects 
taking time to accumulate, so solutions of the perturbed orbit equations would provide a 
reasonable picture of the early development of tidal structures (see review of Struck 1999). 
This, in turn, would be a very efficient means of constraining the satellite orbit and the pri¬ 
mary gravitational potential needed to produce observed tidal features, at least in encounters 
with modest dynamical friction. 

The second is based on the slowly precessing orbits in the inner regions of the softened 
logarithmic potential discussed in Sec. 3.3 and illustrated in Figure 3a. Equation (15) 
tells us how the precession b depends on orbit size p in this case. In particular, for values 
of p < 1, we have &i ~ —1, which is the slow precession. If a disturbance (e.g., tidal) 
elongates initially circular orbits in this region, the slow precession tells us that they will 
stay approximately aligned for some time. This in itself is not remarkable in a potential that 
is approximately solid-body in that region. However, the p-cllipse equations not only tell 
us about this kinematic bar, but they also tell us quantitatively how this structure would 
precess into a spiral at somewhat larger radii. Of course they do not include the effects of 
self-gravity on these kinematic waves. 

A third application concerns resonant orbit-satellite interactions. Tremaine & Wein¬ 
berg (1984), for example, have emphasized the importance of resonant orbits in galaxy bar 
dynamics and in dynamical friction effects on satellites. Frequently, when resonances are 
discussed in galaxy dynamics they are Lindblad resonances. These are local (in radius), and 
usually visualized in the epicyclic approximation as small integer ratios between the epicyclic 
period and the period of the mean (guiding center) orbit. 

The existence of closed orbits and nearly closed orbits described above naturally raises 
the topic of resonant interactions in galaxy disks involving these orbits. One example would 



be resonances between eccentric p-ellipse type orbits and circular orbits with resonant pe¬ 
riods, and radii close to the apsides of the eccentric orbit. Because of the radial range of 
the eccentric orbit, such interactions would not be as localized as the Lindblad resonances. 
For non-Keplerian potentials, Kepler’s third law generally has a (complicated) dependence 
on eccentricity (see Appendix C), which will make the cataloguing of such resonances some¬ 
what difficult. These resonances seem qualitatively similar to those between orbit families 
in barred potentials, so like the latter they may play a role in building bulges, especially at 
early times. 

Period resonances involving the precession of p-ellipse orbits may also be important. 
Figures 6a,b remind us that potentials close to the logarithmic potential with such orbits 
have precession periods of about a few times the orbital period. In these nearly flat rotation 
curve potentials, this means that the precession periods of disk stars are about equal to the 
orbital periods of satellites with orbital radii a few times larger. Higher order resonances 
would occur at larger satellite orbital radii. Thus, especially in the early stages of galaxy 
evolution, this interaction could play a very important role in the evolution of disks and 
satellite orbits. It may also be important in galaxy interactions. The immediate objection is 
that the closed orbit potentials described above are very sparsely distributed on the 5 axis, 
i.e., very rare. A counter-argument is that potentials with orbits that are nearly closed, i.e., 
do not drift much in phase over a duration of several precession periods, span finite intervals 
on the 5 axis near the closed orbit potentials. 

Resonant orbit-precession interactions might also be important between forming planets 
in accretion disks if the disk mass is great enough to modify the Kepler potential, as seen by 
planets on elliptical-like orbits. As described above, there are many closed orbit potentials 
with 5 slightly greater than the Kepler value of 1/2. The precession rate is slow for these 
potentials, but the evolutionary times extend over at least thousands of orbital periods, so 
this is not a problem. 

A similar application may be found for the orbits of stars in dense clusters around 
massive black holes in galaxy nuclei. If the central black hole dominates, but the distributed 
mass of the star cluster is non-negligible, the effective potential will have a rotation curve 
that is slightly flatter than a Keplerian one. This is again in the region of dense closed orbit 
potentials, so precessional resonance interactions with satellite star clusters or other massive 
black holes may be important. As usual with resonant interactions, they likely involve only 
a small minority of stars, and so, specially designed N-body simulations are needed to study 
the phenomenon. 
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6. Summary and Conclusions 

In the sections above, we have shown that curves derived as powers of ellipse functions, 
called p-ellipses, provide simple, yet surprisingly accurate approximations to orbits in a range 
of power-law potentials, including the well-known logarithmic potential, and softened power- 
law potentials. For a given power-law potential the p-cllipse function is nearly as simple as 
an ellipse, but the family of p-cllipses extends continuously across a physically interesting 
range of potentials. 

There are at least two reasons why p-ellipses are good orbital approximates across this 
range of potentials, and are likely to be the best simple, analytic functions to do so. First, 
p-ellipses match the tendency for orbits of a given energy and angular momentum to become 
more circular as the potential changes from Keplcrian to solid-body (and 5 decreases). In 
this sense the p-ellipses adjust well to the appropriate orbital shape in a given potential. 
Secondly, the precession rate obtained by demanding that the p-ellipse satisfy the equation 
of motion in a given potential to first order in eccentricity is the identical to that given 
by Newton’s theorem. Thus, the p-ellipses precess correctly. Moreover, a second order 
approximation yields an eccentricity dependence of the precession rate that is in qualitative 
agreement with the Valluri et al. (2005) semi-analytic results. 

A number of the results obtained above for p-ellipses, like the apsidal precession rates 
in power-law potentials, are not new, and good approximations to individual orbits can be 
obtained with epicycles or numerical integration. However, a family of simple curves like 
the p-cllipses allow us to readily see trends across a range of potentials, so they provide a 
simple, conceptual picture for orbital variations, as discussed in the introduction. 

Moreover, as shown in the previous sections the p-cllipses provide a very powerful tool 
for studying characteristics like the occurrence of (non-circular) closed orbits as a function 
of eccentricity in different potentials. The key feature of the p-ellipses in this regard is a 
relatively simple approximate formulae for precession rates as a function of 6 and e. Newton’s 
theorem is valid in the limit of small eccentricity, and Valluri et al.’s extensions involve 
integrals that must be evaluated numerically. Equations (19) and (20), though approximate, 
offer convenience. (Also see eq. (25).) 

The example of the softened power-law potentials holds out the hope that the p- 
cllipses can provide useful orbit approximations in other non-power-law potentials. Non- 
axisymmetric potentials have not been considered in this paper, but it seems reasonable to 
hope that p-cllipses could provide good approximations to loop orbits in such potentials. 
This issue deserves more study. 

In Section 5.3 I described several examples of how the systematics of p-cllipses might 
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shed light on important astrophysical problems. This conclusion will be even more general 
if p-cllipses prove to be good orbital approximations in more types of potential. 

Orbit theory is more general than celestial mechanics, and the p-cllipse approximations 
should also be relevant to any held than involves orbits in general potentials. Electron orbits 
in general, steady (gradient) electric and magnetic fields are obvious examples. 

In sum, there seems to be room for a great deal more development of the theory and 
application of these simple curves. Their simplicity may allow us to address a number of 
complex issues that would be hard to study directly through the accumulation of numerical 
examples. 
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A. Some Perturbation Analysis Details 

The purpose of this brief appendix is to provide a few more details, and a slight extension 
of the perturbation analysis leading to the results of equation (10). This is the simplest of 
several such calculations in the main text, but is representative of the common procedure. 

That procedure consists of substituting the adopted solution (first equality of eq. (10)) 
into the equation of motion (8), and gathering terms of common order in the factor e cos((l — 
6)0). Then the requirement that these terms cancel up to a given order is used to derive 
relations between the paramenters. In this case, the zeroth order terms yield the condition, 


1 




(1 — 6) 2 e 2 \ 

4 )’ 


(Al) 


for the constant c of equation (8) in terms of the variables b , p, and e. Note that since this 
is an expansion in terms of ecos((l — 6)0), the e 2 term (with no cosine factor) is included 
in the zeroth or constant condition. This factor is not included in the (first-order) results 
given in equation (10). 

Terms of first order in ecos((l — 6)0) yield the condition, 


1(1 - bj 2 - 2 + cp 2 = 0, 


(A2) 
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which, in combination with the previous condition, yield an expression for the b factor. 
The terms of second order in ecos((l — b)(f>) are, 


^ (1 - b) 2 _ 1 e 2 

4 “ 2 ~~ T’ 

where a common factor of 1/p 2 has been dropped, and the result for (1 — 6) from the previous 
condition has been used. With no extra parameter to adjust, these terms do not cancel, but 
do provide an estimate of the p-ellipse accuracy. It is interesting that while all second order 
terms are small at small eccentricity, the term above is smallest at large eccentricity. This 
is not the case for other power-law potentials. 



B. Relations Between Orbital Parameters 

In this appendix I present some relations between orbital parameters, especially those 
between the semi-major axis a, the eccentricity parameter e, and the specific angular mo¬ 
mentum h and specific energy £. The relation between h and the constant c is given in 
equation (4), and equation (19) relates the latter to the semi-latus rectum p, for a first-order 
p-ellipse. Combining these yields, 


h 2 

GM * 


= p 2 ^ = 


9o(e) 


2 ( 1 - 5 ) 


with, g 0 (e) = 


1 


(1 - e) 1/2+<5 + (1 + e) 1/2+<5 

(1 - e 2)!/2+5 


(Bl) 

(B2) 


The second equality in (Bl), h versus a, is derived from a generalization of equation (11) for 
the case <5^0. 

The specific energy is defined as, 


and, £ 



= i(u 2 -GM*ln(r)), 


for <5^0, 

(B3) 

for 8 = 0, 

(B4) 


where v is the total relative velocity of satellite. The radial velocity, v r , can be derived 
in terms of the azimuthal velocity by differentiating equation (18) with respect to time. 
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Then, the azimuthal velocity can be eliminated using the relation h = rv After some 
simplification, for the case with 5 ^ 0, we get, 


with, g e (e,5) 


g 2 0 5 ( 1 + e) 


1+2 8 


c GM * ( 

^ 2 a 26 


1 - 4(1 + e)(‘ 5_1 )d+ 2<5 ) 
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(B5) 
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For the case <5 = 1/2, equations (Al), (A2), (A5), and (A6) reduce to the usual Keplerian 
elliptical orbit relations: £ = —GM/2a and h 2 = GMa( 1 — e 2 ). 

For the case 0 < <5 <C 1/2 (i.e., nearly flat, but slightly decreasing rotation curve), we 
have, 


h 2 


9o - 


= GM*a 2 ( 1 - 
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g £ ~l + e--~ 
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(B7) 
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where e Q is the “classical” eccentricity of equation (23). In contrast to the full relations, 
these limiting cases can be easily inverted to derive a, e as functions of h, £. 


C. Time-Azimuth Relations and Kepler’s Law for p-ellipses in Power-Law 

Potentials 

Relations between time and azimuth for p-cllipse orbits in power-law potentials are 
not treated in the main text. Since these relations are likely to be very important for any 
application involving p-cllipses we derive them here. Specifically, we can derive an invertible 
analytic expression for the logarithmic potential, and a implicit relation for potentials with 
small values of <5. To begin, the equation for the conserved specific angular momentum gives 
the azimuthal frequency, d(j)/dt = h/r 2 . Next, we substitute the p-ellipse solution for r to 
get, 


d(f) 

dt 


■4 [1 + e cos ((1 - b)(j))) 1+25 . 
jr 


(Cl) 


In the case of the logarithmic potential, <5 = 0, this can be integrated to obtain, 
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1 — e 
1 + e 


tan I (1 ~ 6) ^ - tan < (1 “ ^ 


(C2) 


where, on the right hand side, we assume that there is an initial azimuth, 0 O , on the left we 
assume that the initial time is t — 0. 

By setting (1 — b)(f> = 2ix and 0 O = 0 we get an expression for the orbital period, T, in 
the logarithmic potential, 


T 


p 2 2vr 

h( 1 - b) yj\ - e 2 ’ 


(C3) 


Using equations (10) for b , (11) for p in terms of the semi-major axis a, and (Bl), 
(B2) for h (with <5 = 0), we derive Kepler’s third law for first order p-cllipse orbits in the 
logarithmic potential, 


rp2 


2vr 2 

GM* 


n 2 


(a/1 + e + a/1 - e) J 


a . 


(C4) 


The eccentricity dependence of the expression in square brackets is less than 10% different 
than the e = 0 value when e < 0.7, and so is quite weak. 

For power-law potentials near the logarithmic, |h| -C 1, the procedure above gives the 
following integral equation, 


ht 


dcf) 1 

[1 + ecos ((1 — b)(f>)} 1+25 1 — b 


[1 — 2 Se cos ((1 — b)4>)\ 
[1 + e cos ((1 — b)(p)] 


d(( 1 - b)<j >), 


(C5) 


where the last equality gives a first order approximation in 5. This equation can be further 
reduced to, 


ht 
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1-26 

1-6 


I 0 -26(</>- <j> 0 ), with, I a = 


d(j) 


Ko [1 +ecos ((1-6)0)] 


(C6) 


The integral I 0 is the same one that appears in the case of the logarithmic potential, and 
whose solution is give in equation (C2). The above equation cannot be solved explicitly for 
0. However, an explicit form of Kepler’s third law can be derived; it is, 



2 


rj ~>2 


2vr 2 

'1- 26- JVl-e 2 " 
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2(1 — e 2 ) 5 

GM* 

L vi - J 


(jl + e)3+ 5 + (1 - e)^ +<5 ) 


a 2{ - l+5 \ 


(C7) 


For small positive values of 6, this expression gives an e dependence that is not quite as weak 
as that of the logarithmic potential. 
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Fig. 1.— Sample x-y orbit in the logarithmic potential with energy/angular momentum 
constant c = 1.2 and eccentricity e = 0.3. The solid line is a numerical integration of the 
orbit, dashed line is a first-order p-ellipse approximation. See text for details. Note the end 
points of both in the lower right quadrant. 
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Fig. 2.— Inverse radius variable u versus azimuthal angle 0 for a sample orbit in the 
logarithmic potential with energy/angular momentum constant c = 1.2 and eccentricity 
e = 0.3. The solid line is a numerical integration of the orbit, dashed line is the first-order 
p-ellipse approximation, and dotted line is the second-order approximation. 
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Fig. 4.— Four sample orbits in several power-law potentials with exponent 6 values of 
S = —0.25, 0.25, 0.45, 0.75. The p-ellipse eccentricity e = 0.2 and orbit size p — 1 in all cases. 
The solid line is a numerical integration of the orbit, dashed line is the first-order p-ellipse 
approximation. 
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Fig. 5.— Two sample orbits in power-law potentials with steeply rising rotation curves, i.e., 
exponent 5 values of 5 = —1.1,—1.5. The p-cllipse eccentricity e = 0.2 and orbit size p = 1 
as in the previous figure. The solid line is a numerical integration of the orbit, dashed line 
is the first-order p-cllipse approximation. 
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Fig. 6.— Four sample closed orbits in power-law potentials with exponent 6 values of 5 = 
— 1/8,1/9, 7/32, 0.44875. The p-ellipse eccentricity e = 0.2 and orbit size p = 1 in all cases. 
The solid line is a numerical integration of the orbit, dashed line is the first-order p-cllipse 
approximation. In the lower right panel only the numerical curve is shown for clarity, though 
the approximation agrees as well as in the other cases. 
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Fig. 7.— Inverse radius variable u versus azimuthal angle 0 for a sample orbit in the 
5 = 1/9 potential with energy/angular momentum constant c = 1.0 and eccentricity e = 0.95. 
The solid line is a numerical integration of the orbit, dashed line is the first-order p-cllipse 
approximation, and dotted line is the second-order approximation. 
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Fig. 8.— Sample x-y orbit in the 5 = 1/9 potential with energy/angular momentum constant 
c = 1.0 and eccentricity e = 0.95. The solid curve is a numerical integration of the orbit, 
dashed curve is the first-order p-ellipse approximation, and dotted curve is a spirographic 
approximation as described in the text. 
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Table 1: Examples of Potentials with Closed Orbits of Low Eccentricity. 
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See text for definitions of b and <5 



